Basics of structural equation modeling
04.01.2024
Dataset
The data for following exercises stem from the publication Grassland ecosystem recovery after soil disturbance depends on nutrient supply rate (Seabloom, Borer, and Tilman (2020)) and are publicly available at Dryad. The data were obtained during the long-term field experiment Cedar Creek LTER and target the effects of human disturbances on grassland ecosystem functioning and biodiversity.
Data collection
Following description of the data collection process is taken from Seabloom, Borer, and Tilman (2020).
Treatment
The disturbance treatment was replicated in the three old-fields (A, B and C) in a completely randomised block design (two treatments in each of three fields for a total of 6 35 × 55 m large plots). In 1982, in each of the fields, one of these two 35 × 55 m areas was selected to be disturbed with a disk harrow. After, the soil was hand raked to smooth the soil and remove any remaining vegetation, so that subsequent colonisation was solely from seeds or small rhizome fragments. Within each of the 6 large plots, the 54 small plots were arrayed in 6 × 9 grid with 1 m buffers between each plot. Aluminium flashing was buried to depth of 30 cm around each plot to prevent horizontal movement of nutrients and spreading of plants through vegetative growth.
The nutrient treatments were replicated six times in a completely randomised design in each of the 35 × 55 m plots (54 4 × 4 m small plots) yielding 324 (6 x 54) plots. The Seabloom analysis focuses on two nutrient treatments:
- Control (no nutrients; Treatment I) and
- Other Nutrients and 9.5 g of N (Treatment F)
Sampling and analysis
At peak biomass (mid-July to late August), all aboveground biomass was clipped in a 3 m by 10 cm strip (0.3 m2) in each plot. Note that there were 4 years when the disturbed plots were not sampled or only sampled in a single field. The biomass was sorted into dead, previous year’s growth (litter) and current year’s growth (live biomass). Live biomass was sorted to species, dried and weighed. We estimated total aboveground biomass as the summed biomass of all non-woody species in each 0.3 m2 sample, converted to g/m2.
Species richness is the number of species in each 0.3 m2 sample. We quantified plant diversity as the Effective Number of Species based on the Probability of Interspecific Encounter (ENSPIE), a measure of diversity that is more robust to the effects of sampling scale and less sensitive to the presence of rare species than species richness. ENSPIE is equivalent to the Inverse Simpson’s index of diversity (\(1 / \sum_{i=1}^{S} p_i^2\) where \(S\) is the total number of species and \(p_i\) is the proportion of the community biomass represented by species \(i\)).
Load data
seabloom <- read.table(here("2_Modeling/Data_preparation/seabloom-2020-ele-dryad-data/cdr-e001-e002-output-data.csv"),
sep = ",", header = TRUE)Explore data
dim(seabloom)## [1] 5040 16str(seabloom)## 'data.frame': 5040 obs. of 16 variables:
## $ exp : int 1 1 1 1 1 1 1 1 1 1 ...
## $ field : chr "A" "A" "A" "A" ...
## $ plot : int 1 1 1 1 1 1 1 1 1 1 ...
## $ disk : int 0 0 0 0 0 0 0 0 0 0 ...
## $ yr.plowed : int 1968 1968 1968 1968 1968 1968 1968 1968 1968 1968 ...
## $ ntrt : int 9 9 9 9 9 9 9 9 9 9 ...
## $ nadd : num 0 0 0 0 0 0 0 0 0 0 ...
## $ other.add : int 0 0 0 0 0 0 0 0 0 0 ...
## $ year : int 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 ...
## $ dur : int 1 2 3 4 5 6 7 8 9 10 ...
## $ precip.mm : num 879 858 876 994 859 ...
## $ precip.gs : num 374 551 527 503 512 ...
## $ mass.above: num 61.3 135.9 216.9 302.8 586.9 ...
## $ rich : int 13 15 14 14 16 14 8 7 9 13 ...
## $ even : num 0.463 0.156 0.204 0.14 0.269 ...
## $ ens.pie : num 6.02 2.34 2.85 1.96 4.31 ...Exploring the data reveals 16 variables with each 5040 data points:
exp: treatments in split-plot design: 1 = disturbance (Control or Disked, 35 × 55 m plots) and 2 = nutrient addition (9 levels, 4 × 4 m plots)field: three experimental fields A, B and Cplot: 54 plots within fieldsdisk: disking treatment (0 = intact at start of experiment, 1 = disked at start of experiment)yr.plowed: last year field was plowed for agriculture (A: 1968, B: 1957 and C: 1934)ntrt: nine levels representing different combinations of nitrogen (0 to 27.2 g N year-1 added as NH4NO3) and other nutrients (20 g m−2 year−1 P205; 20 g m−2 year−1 K20; 40 g m−2 year−1 CaCO3; 30.0 g m−2 year−1 MgSO4; 18 μg m−2 year−1 CuSO4; 37.7 μg m−2 year−1 ZnSO4; 15.3 μg m−2 year−1 CoCO2; 322 μg m−2 year−1 MnCl2 and 15.1 μg m−2 year−1 NaMoO4; details see Table S1 in publication). Nutrients were applied twice per year in mid-May and mid-June.nadd: nitrogen additon rate (g/m2/yr)other.add: other nutrient treatment (0 = control, 1 = other nutrients added)year: sampling yeardur: duration of experimentprecip.mm: annual precipitation (mm)precip.gs: growing season precipitation (mm)mass.above: aboveground biomass (g/m2)rich: species richness (species/0.3 m2)even: Simpson’s evennessens.pie: effective number of species (= probability of interspecific encounter, equivalent to inverse Simpson’s diversity)
Factors as factors
As ntrt is a coding scheme for different nutrients
added, it is crucial to treat it as such rather than as integers.
seabloom$ntrt <- as.factor(seabloom$ntrt)
str(seabloom$ntrt)## Factor w/ 8 levels "2","3","4","5",..: 8 8 8 8 8 8 8 8 8 8 ...Exercise: what do you think is wrong with keeping
ntrt as integers?
Subset dataset
To simplify the analysis, we will take the average across all the years in the dataset for our analysis.
seabloom <- seabloom %>% group_by(exp, field, plot, disk, yr.plowed, ntrt, nadd, other.add) %>% summarise(across(mass.above:ens.pie, mean))## `summarise()` has grouped output by 'exp', 'field', 'plot', 'disk',
## 'yr.plowed', 'ntrt', 'nadd'. You can override using the `.groups` argument.Overview
The pairs function yields an overview over the numerical
data that we will use for the following exercises.
ggpairs(ungroup(seabloom[, c(7, 9:12)]), upper=NULL, lower = list(continuous = wrap("points", alpha = 0.05))) + theme_bw() + theme(panel.grid=element_blank())
Metamodel
A metamodel summarizes the concept behind a model and links it to theory. Here, the metamodel is visualized as a directed acyclic graph (DAG) which reads as: productivity (biomass) is directly influenced on the one hand by the environment (nutrients, disturbance and precipitation) and on the other hand by biodiversity (richness and evenness). Also some elements of the environment influence biodiversity and thus, have an additional indirect effect on productivity via biodiversity.
Linear model
First, implement the metamodel into a linear model (LM). For this, three models are necessary: one that accounts for the direct- and two for the indirect effects.
Direct effects
If a linear model is visualized in a DAG, it becomes apparent that there is a set of permitted, but unanalyzed correlations among the predictors. These correlations have a huge influence on the coefficients between the predictors \(x_{1\ldots i}\) (here, nutrients, precipitation, richness and evenness) and the response \(y\) (here, biomass). Despite their importance, it is impossible to include any information about the reason of the correlations between the predictors. Further, they make it nearly impossible to create a proper causal model, since there are many possible causal relations that can create a set of “unanalyzed associations” which hinders interpretations (James B. Grace 2021).
Exercise: write the linear model in the DAG above in R syntax.
# lm.dir <- lm(mass.above ~ nadd + rich + even, data = seabloom)
# summary(lm.dir)Indirect effects
To account for the indirect effects, two additional LMs are necessary: one with richness and one with evenness as response.
lm.rich <- lm(rich ~ nadd, data = seabloom)
summary(lm.rich)##
## Call:
## lm(formula = rich ~ nadd, data = seabloom)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.7062 -1.5628 -0.4956 0.9744 8.9675
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.42278 0.18432 45.70 <2e-16 ***
## nadd -0.19025 0.01528 -12.45 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.298 on 286 degrees of freedom
## Multiple R-squared: 0.3516, Adjusted R-squared: 0.3494
## F-statistic: 155.1 on 1 and 286 DF, p-value: < 2.2e-16lm.even <- lm(even ~ nadd, data = seabloom)
summary(lm.even)##
## Call:
## lm(formula = even ~ nadd, data = seabloom)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.203129 -0.051110 -0.004937 0.046817 0.278873
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3564682 0.0063701 55.96 <2e-16 ***
## nadd 0.0057968 0.0005279 10.98 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07941 on 286 degrees of freedom
## Multiple R-squared: 0.2966, Adjusted R-squared: 0.2941
## F-statistic: 120.6 on 1 and 286 DF, p-value: < 2.2e-16Conclusion
The direct effect model showed that biomass was statistically
significantly positively related to nutrient input (the nitrogen additon
rate nadd) and species richness (rich), but
negatively to evenness (even). Additionally, nutrient input
had a statistically significantly negative influence on species richness
but a positive effect on evenness.
Structural equation modeling
Structural equations aspire to represent cause-effect relationships and thus, they can be used to represent scientific, causal hypotheses. A key feature of structural equation models is the ability to investigate the networks of connections among system components (James B. Grace et al. 2012).
To evaluate the SEMs, we will use the package lavaan that relies on
the computation of covariance matrices to fit the structural equations
(this approach is known as ‘global estimation’). It comes with full
support for categorical data (any mixture of binary, ordinal and
continuous observed variables), can handle both latent and composite
variables, performs mediation analysis and calculates the magnitude of
indirect effects (Rosseel 2012, 2021).
library("lavaan")## This is lavaan 0.6-17
## lavaan is FREE software! Please report any bugs.Collinearity
Inspecting our set of variables shows that correlations between all variables are rather weak except the precipitation throughout the year and those throughout the field season.
ggpairs(ungroup(seabloom[, c(7, 9:12)]), lower = list(continuous = wrap("points", alpha = 0.05))) + theme_bw() + theme(panel.grid=element_blank())
Variables that are highly correlated with \(r > 0.85\) become redundant. Then, one of them can be dropped or they can be modeled summarized as a latent variable (James B. Grace 2006).
An example SEM
This toy model shall illustrate the logic of SEM. It contains the same variables and aims at evaluating the same metamodel as the LM before to ease comparability between the two methods. A huge benefit of SEM is that variables can appear as both predictors and responses what allows the evaluation of direct and indirect effects in one go (as shown in the directed acyclic graph (DAG) below). In SEM jargon, predictors (nodes that have no arrows pointing at them) are called exogenous-, while responses (nodes that have arrows pointing at them) endogeneous variables.
The table below summarizes the operators of lavaan
syntax. For example, an arrow in a DAG is represented by a tilde.
| Formula type | Operator | Meaning | Example |
|---|---|---|---|
| regression | ~ |
is regressed on | y ~ x |
| correlation | ~~ |
correlate errors for | y1 ~~ y2 |
| latent variable | =~ |
set reflective indicators | Height =~ y1 + y2 + y3 |
| composite variable | <~ |
set formative indicators | Comp1 <~ 1*x1 + x2 + x3 |
| mean/intercept | ~ 1 |
estimate mean for y without
xs |
y ~ 1 |
| label parameter | * |
name coefficients | y ~ b1*x1 + b2*x2 |
| define quantity | := |
define quantity | TotalEffect := b1*b3 + b2 |
| define thresholds | | |
for ordered categorial variable | u | t1 + t2 |
Exercise: translate the model in the DAG above into
lavaan syntax with the help of this table.
simple <-
"mass.above ~ nadd + rich + even + disk
rich ~ nadd
even ~ nadd
"Normal distribution
lavaan uses a \(\chi^2\) test to compare the estimated-
to the observed covariance matrix to compute the goodness of fit for the
SE model under the assumption that all observations are independent and
all variables follow a (multivariate) normal distribution (James B. Grace
2006). Note, that these distributional assumptions only apply
to endogenous variables, whereas the distribution of exogenous variables
has no bearing on assumptions (James B. Grace 2021).
We use a graphical method to assess the fit of the endogeneous variables to a normal distribution, the quantile-quantile plots (Q-Q plots). Hereby, the quantiles of the data are compared to those of a theoretical distribution (i.e., the normal distribution). If the data would be normally distributed, the points would match one to one and thus, align diagonally. In the Q-Q plot, this match is indicated by the black line.
The package MVN
allows to plot several Q-Q plots at once and also offers several tests
to multivariate
normal distribution (MVN). Here, we employ the Henze-Zirkler’s test
that has been recommended as a formal test of MVN (Mecklin and Mundfrom
2005).
library("MVN")
mvn(data = seabloom[, c(9:11)], mvnTest = "hz", univariatePlot = "qqplot")## $multivariateNormality
## Test HZ p value MVN
## 1 Henze-Zirkler 5.06108 0 NO
##
## $univariateNormality
## Test Variable Statistic p value Normality
## 1 Anderson-Darling mass.above 5.9190 <0.001 NO
## 2 Anderson-Darling rich 3.5448 <0.001 NO
## 3 Anderson-Darling even 1.9401 1e-04 NO
##
## $Descriptives
## n Mean Std.Dev Median Min Max
## mass.above 288 272.2597298 95.7125006 260.9607974 139.2747737 662.4553000
## rich 288 6.8651316 2.8486571 6.4473684 2.0000000 17.2000000
## even 288 0.4039294 0.0945156 0.3904727 0.2274666 0.7930136
## 25th 75th Skew Kurtosis
## mass.above 200.5621772 318.6579392 1.3912202 2.6381413
## rich 4.7236842 8.3684211 0.8675777 0.5760362
## even 0.3445196 0.4512502 0.8754018 1.5768164
Alternatively, also histograms overlaid with a normal distribution of the same mean and standard deviation as the data allows to gain insight into the distribution. First, let’s define a function that plots the histogram and the expected density curve (it expects a numerical vector and a character string as input):
histWithDensity <- function(variable, name){
hist(variable, prob = TRUE, main = "", xlab = name)
x <- seq(min(variable), max(variable), length = 400)
y <- dnorm(x, mean = mean(variable), sd = sd(variable))
lines(x, y, col = "red", lwd = 2)
}Then, we can apply this function to our three endogenous variables to plot the histograms with the expected density curves under a normal distribution:
par(mfrow = c(1, 3))
histWithDensity(seabloom$mass.above, "mass.above")
histWithDensity(seabloom$rich, "rich")
histWithDensity(seabloom$even, "even")
par(mfrow = c(1, 1))Exercise: would you infer that the endogenous variables meet the assumption of being (multivariate) normally distributed from the results of the the Q-Q plots, the Henze-Zirkler test and the histograms? And why (not)?
Data transformation
Transformation of the three endogenous variables shall make them more “normal”.
seabloom$log.mass.above <- sqrt(seabloom$mass.above)
seabloom$log.even <- sqrt(seabloom$even)
seabloom$log.rich <- sqrt(seabloom$rich)
mvn(data = seabloom[, c(13:15)], mvnTest = "hz", univariatePlot = "qqplot")## $multivariateNormality
## Test HZ p value MVN
## 1 Henze-Zirkler 3.205733 2.220446e-16 NO
##
## $univariateNormality
## Test Variable Statistic p value Normality
## 1 Anderson-Darling log.mass.above 2.7243 <0.001 NO
## 2 Anderson-Darling log.even 0.7415 0.0528 YES
## 3 Anderson-Darling log.rich 0.8080 0.0362 NO
##
## $Descriptives
## n Mean Std.Dev Median Min Max
## log.mass.above 288 16.2747569 2.7235587 16.1542806 11.8014734 25.7382070
## log.even 288 0.6314047 0.0726349 0.6248782 0.4769346 0.8905131
## log.rich 288 2.5661429 0.5301115 2.5391616 1.4142136 4.1472883
## 25th 75th Skew Kurtosis
## log.mass.above 14.1619976 17.8509923 0.8659716 0.8744144
## log.even 0.5869577 0.6717516 0.4606483 0.5892468
## log.rich 2.1733975 2.8928223 0.3724584 -0.2250292par(mfrow = c(1, 3))
histWithDensity(seabloom$log.mass.above, "log(mass.above)")
histWithDensity(seabloom$log.rich, "log(rich)")
histWithDensity(seabloom$log.even, "log(even)")
par(mfrow = c(1, 1))Fit the model
Now, let’s fit the model with lavaan’s sem
function. As the data clearly deviates from a MVN, we use the
MLM estimator that provides standard errors and a \(\chi^2\) test statistic robust to
non-normality. Hereby, the Satorra-Bentler correction is used to correct
the value of the ML-based \(\chi^2\)
test statistic by an amount that reflects the degree of kurtosis (Rosseel 2012).1
fit.simple <- sem(simple, data = seabloom, estimator = "MLM")## Warning in lav_data_full(data = data, group = group, cluster = cluster, :
## lavaan WARNING: some observed variances are (at least) a factor 1000 times
## larger than others; use varTable(fit) to investigateOups, the algorithm converges with a warning. Kindly, it informs us how to fix this.
Exercise: let’s obey the software and execute the code from the hint.
# varTable(fit.simple)This reveals an enormous difference in magnitude between the variance
of biomass (mass.above) and the other variables.
Rescale variables
To remove this difference in magnitude between the variables, we
divide mass.above and precip.mm by 100 what
changes the unit from g/m2 to 10 mg/m2 and from mm
to 0.1 m respectively. The boxplots show that the range of the variables
is now more similar.
seabloom$mass.above <- seabloom$mass.above / 100
boxplot(seabloom[, c("mass.above", "rich", "even", "nadd", "disk")], las = 2)
Then, run sem again with the rescaled
biomass variable:
fit.simple <- sem(simple, data = seabloom, estimator = "MLM")
summary(fit.simple, fit.measures = TRUE)## lavaan 0.6.17 ended normally after 1 iteration
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 9
##
## Number of observations 288
##
## Model Test User Model:
## Standard Scaled
## Test Statistic 120.437 118.894
## Degrees of freedom 3 3
## P-value (Chi-square) 0.000 0.000
## Scaling correction factor 1.013
## Satorra-Bentler correction
##
## Model Test Baseline Model:
##
## Test statistic 592.561 765.397
## Degrees of freedom 9 9
## P-value 0.000 0.000
## Scaling correction factor 0.774
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.799 0.847
## Tucker-Lewis Index (TLI) 0.396 0.540
##
## Robust Comparative Fit Index (CFI) 0.800
## Robust Tucker-Lewis Index (TLI) 0.399
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -597.880 -597.880
## Loglikelihood unrestricted model (H1) -537.662 -537.662
##
## Akaike (AIC) 1213.761 1213.761
## Bayesian (BIC) 1246.728 1246.728
## Sample-size adjusted Bayesian (SABIC) 1218.187 1218.187
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.369 0.366
## 90 Percent confidence interval - lower 0.314 0.312
## 90 Percent confidence interval - upper 0.426 0.424
## P-value H_0: RMSEA <= 0.050 0.000 0.000
## P-value H_0: RMSEA >= 0.080 1.000 1.000
##
## Robust RMSEA 0.369
## 90 Percent confidence interval - lower 0.314
## 90 Percent confidence interval - upper 0.427
## P-value H_0: Robust RMSEA <= 0.050 0.000
## P-value H_0: Robust RMSEA >= 0.080 1.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.110 0.110
##
## Parameter Estimates:
##
## Standard errors Robust.sem
## Information Expected
## Information saturated (h1) model Structured
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## mass.above ~
## nadd 0.096 0.008 11.391 0.000
## rich 0.106 0.016 6.711 0.000
## even -1.378 0.555 -2.486 0.013
## disk 0.446 0.073 6.066 0.000
## rich ~
## nadd -0.190 0.014 -13.573 0.000
## even ~
## nadd 0.006 0.001 8.503 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .mass.above 0.389 0.049 7.933 0.000
## .rich 5.243 0.536 9.786 0.000
## .even 0.006 0.001 10.754 0.000Goodness of fit
This time, the model converged, however, with poor fit:
- The ratio of the test statistic and the degrees of freedom should be smaller than 2. Here, the ratio is 120.44 / 3 = 40.15, which indicates that the model is quite far away from a decent fit (James B. Grace 2021).
- The \(p\)-value, which represents the probability of the data given our model (or no significant discrepancy between model and data), should be larger than 0.05 (James B. Grace 2021).
- The comparative fit index (CFI) ranges between zero and one, whereas a value \(> .95\) is considered as a good fit indicator (Hu and Bentler 1999). Our robust CFI of 0.8 is fairly below this threshold. In comparison to the \(p\)-value, the CFI has the advantage to be independent of sample size.
Modification indices
To improve the model fit, we look for missing paths via the modification indices. They indicate an estimated drop in the model \(\chi^2\) resulting from freeing fixed parameters via including a missing path. 3.84 is considered as the critical threshold, the “single-degree-of-freedom \(\chi^2\) criterion” (James B. Grace 2021).
modindices(fit.simple, minimum.value = 3.84)## lhs op rhs mi epc sepc.lv sepc.all sepc.nox
## 15 rich ~~ even 96.523 -0.105 -0.105 -0.579 -0.579
## 16 rich ~ mass.above 32.400 3.096 3.096 1.017 1.017
## 17 rich ~ even 96.523 -16.751 -16.751 -0.556 -0.556
## 19 even ~ mass.above 55.110 -0.105 -0.105 -1.042 -1.042
## 20 even ~ rich 96.523 -0.020 -0.020 -0.603 -0.603In the column mi (for modification index) we look for
high values. Note, however, that the modification indices are uninformed
suggestions and further adaptations of the model based on their
information needs to be based on theory.
In this example, the modification indices indicate–amongst other–a missing relation between richness and evenness. Including this relation into the model would expectedly improve its fit by a change in the \(\chi^2\) by 96.523.
This path is necessary as richness and evenness are computationally
related to each other (they are not independent quantities). Thus, let’s
include a correlation between rich and even
(Note: in SEM, a correlation between two variables points to an
omitted common cause/variable that drives this correlation). With the
function update() it is possible to directly incorporate
the missing path into the specified model without rewriting it from
scratch:
fit.simple.up <- update(fit.simple, add = "rich ~~ even")
summary(fit.simple.up, fit.measures = TRUE, rsq = TRUE)## lavaan 0.6.17 ended normally after 29 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 10
##
## Number of observations 288
##
## Model Test User Model:
## Standard Scaled
## Test Statistic 2.878 2.892
## Degrees of freedom 2 2
## P-value (Chi-square) 0.237 0.235
## Scaling correction factor 0.995
## Satorra-Bentler correction
##
## Model Test Baseline Model:
##
## Test statistic 592.561 765.397
## Degrees of freedom 9 9
## P-value 0.000 0.000
## Scaling correction factor 0.774
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.998 0.999
## Tucker-Lewis Index (TLI) 0.993 0.995
##
## Robust Comparative Fit Index (CFI) 0.998
## Robust Tucker-Lewis Index (TLI) 0.993
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -539.101 -539.101
## Loglikelihood unrestricted model (H1) -537.662 -537.662
##
## Akaike (AIC) 1098.202 1098.202
## Bayesian (BIC) 1134.831 1134.831
## Sample-size adjusted Bayesian (SABIC) 1103.120 1103.120
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.039 0.039
## 90 Percent confidence interval - lower 0.000 0.000
## 90 Percent confidence interval - upper 0.130 0.130
## P-value H_0: RMSEA <= 0.050 0.457 0.454
## P-value H_0: RMSEA >= 0.080 0.302 0.305
##
## Robust RMSEA 0.039
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.130
## P-value H_0: Robust RMSEA <= 0.050 0.456
## P-value H_0: Robust RMSEA >= 0.080 0.302
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.021 0.021
##
## Parameter Estimates:
##
## Standard errors Robust.sem
## Information Expected
## Information saturated (h1) model Structured
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## mass.above ~
## nadd 0.096 0.007 13.956 0.000
## rich 0.106 0.019 5.655 0.000
## even -1.378 0.665 -2.072 0.038
## disk 0.446 0.073 6.066 0.000
## rich ~
## nadd -0.190 0.014 -13.573 0.000
## even ~
## nadd 0.006 0.001 8.503 0.000
##
## Covariances:
## Estimate Std.Err z-value P(>|z|)
## .rich ~~
## .even -0.105 0.011 -9.600 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .mass.above 0.389 0.049 7.933 0.000
## .rich 5.243 0.536 9.786 0.000
## .even 0.006 0.001 10.754 0.000
##
## R-Square:
## Estimate
## mass.above 0.569
## rich 0.352
## even 0.297Now, the model has a decent fit with a ratio of test statistic and degrees of freedom smaller than two (i.e. 1.44) and a \(p\)-value of 0.24. Further, also the robust CFI is now 1.
Another look at the modification indices shows that only one modification would yield an improvement, but based on the experimental design AGB cannot influence the disk treatment which was randomly assigned. All other modifications would yield only smallest improvements to the model fit, so that we can ignore them confidently:
modindices(fit.simple.up, minimum.value = 0.01)## lhs op rhs mi epc sepc.lv sepc.all sepc.nox
## 16 rich ~ mass.above 2.779 0.823 0.823 0.275 0.275
## 18 rich ~ disk 2.779 0.367 0.367 0.064 0.129
## 19 even ~ mass.above 0.529 0.012 0.012 0.125 0.125
## 21 even ~ disk 0.529 0.006 0.006 0.029 0.059
## 26 disk ~ mass.above 6.435 0.472 0.472 0.898 0.898
## 27 disk ~ rich 0.965 0.008 0.008 0.046 0.046
## 28 disk ~ even 0.040 -0.052 -0.052 -0.010 -0.010Model comparison
To evaluate whether the simple or the updated model perform better,
we can calculate the Akaike
information criterion (AIC) and compute a significance test based on
the \(\chi^2\) test statistic with the
anova function.
anova(fit.simple, fit.simple.up)##
## Scaled Chi-Squared Difference Test (method = "satorra.bentler.2001")
##
## lavaan NOTE:
## The "Chisq" column contains standard test statistics, not the
## robust test that should be reported per model. A robust difference
## test is a function of two standard (not robust) statistics.
##
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## fit.simple.up 2 1098.2 1134.8 2.8779
## fit.simple 3 1213.8 1246.7 120.4370 112.08 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1From the ANOVA table we can see the difference in AIC = 115.56 is
clearly in favor for the updated simple.up model
(Note: a difference in the AIC of two is considered as
informative (Burnham and
Anderson 2002)).
Further, the difference in the \(\chi^2\) is 112.08–far beyond the threshold value of 3.84 that is necessary to detect a statistically significant difference on a confidence level of \(\alpha = 0.05\).
Standardized coefficients
Using analysis of covariances allows for estimation of both unstandardized (raw) and standardized coefficients. While the analysis is based on covariances for estimating unstandardized coefficients, it is based on correlations for estimating standardized coefficients. The computational relation between correlations and covariances is:
\(r_{xy} = \frac{cov_{xy}}{SD_x \times SD_y}\)
… where \(r_{xy}\) are the correlations, \(cov_{xy}\) the covariances and \(SD_x\) and \(SD_y\) the respective standard deviations (the square roots of the respective variance \(= \sqrt{var_x}\) and \(= \sqrt{var_y}\)) (James B. Grace 2021).
Significance testing of standardized coefficients violates statistical principles. Significance tests are hence based on unstandardized (raw) coefficients and standardized coefficients are used for interpretation only (James B. Grace 2021).
For raw coefficients, the predicted effects are in raw units and have a pretty straight-forward interpretation: as in our example, e.g., the predicted change in 10 mg/m2 above-ground biomass associated with a certain change in added nutrients. For the standardized coefficients, however, the interpretation is more complex: they express the predicted changes in terms of standard deviation units. Thus, they are only interpretable within a sample which hinders generalization. Still, they ease comparison as they are in the same units across various pathways (James B. Grace 2021).
There are different standardizations available in
lavaan, whereof std.all is based on the
standard deviation:
standardizedsolution(fit.simple.up, type = "std.all")## lhs op rhs est.std se z pvalue ci.lower ci.upper
## 1 mass.above ~ nadd 0.897 0.036 24.781 0.000 0.826 0.967
## 2 mass.above ~ rich 0.316 0.050 6.343 0.000 0.219 0.414
## 3 mass.above ~ even -0.137 0.066 -2.072 0.038 -0.266 -0.007
## 4 mass.above ~ disk 0.235 0.034 6.908 0.000 0.168 0.301
## 5 rich ~ nadd -0.593 0.029 -20.543 0.000 -0.650 -0.536
## 6 even ~ nadd 0.545 0.045 12.116 0.000 0.456 0.633
## 7 rich ~~ even -0.579 0.033 -17.449 0.000 -0.644 -0.514
## 8 mass.above ~~ mass.above 0.431 0.034 12.644 0.000 0.364 0.497
## 9 rich ~~ rich 0.648 0.034 18.939 0.000 0.581 0.715
## 10 even ~~ even 0.703 0.049 14.369 0.000 0.607 0.799
## 11 nadd ~~ nadd 1.000 0.000 NA NA 1.000 1.000
## 12 nadd ~~ disk 0.000 0.000 NA NA 0.000 0.000
## 13 disk ~~ disk 1.000 0.000 NA NA 1.000 1.000Derived quantities
To estimate derived quantities and their standard errors, we can use
the := operator in lavaan. The below code shows how to
derive the direct, indirect and total effect of nutrients on above
ground biomass by calculating compound paths and by summing up direct
and indirect effects.
derived <-
"mass.above ~ b1 * nadd + b2 * rich + b3 * even + disk
rich ~ b4 * nadd
even ~ b5 * nadd
dir.nut.effect := b1
indir.nut.effect := b2 * b4 + b3 * b5
tot.nut.effect := b1 + b2 * b4 + b3 * b5
"
fit.derived <- sem(derived, data = seabloom, estimator = "MLM")
summary(fit.derived, rsq = TRUE)## lavaan 0.6.17 ended normally after 1 iteration
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 9
##
## Number of observations 288
##
## Model Test User Model:
## Standard Scaled
## Test Statistic 120.437 118.894
## Degrees of freedom 3 3
## P-value (Chi-square) 0.000 0.000
## Scaling correction factor 1.013
## Satorra-Bentler correction
##
## Parameter Estimates:
##
## Standard errors Robust.sem
## Information Expected
## Information saturated (h1) model Structured
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## mass.above ~
## nadd (b1) 0.096 0.008 11.391 0.000
## rich (b2) 0.106 0.016 6.711 0.000
## even (b3) -1.378 0.555 -2.486 0.013
## disk 0.446 0.073 6.066 0.000
## rich ~
## nadd (b4) -0.190 0.014 -13.573 0.000
## even ~
## nadd (b5) 0.006 0.001 8.503 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .mass.above 0.389 0.049 7.933 0.000
## .rich 5.243 0.536 9.786 0.000
## .even 0.006 0.001 10.754 0.000
##
## R-Square:
## Estimate
## mass.above 0.554
## rich 0.352
## even 0.297
##
## Defined Parameters:
## Estimate Std.Err z-value P(>|z|)
## dir.nut.effect 0.096 0.008 11.391 0.000
## indir.nut.ffct -0.028 0.006 -4.770 0.000
## tot.nut.effect 0.068 0.006 12.214 0.000Saturated model
In a saturated model all possible paths are specified and as a result, there are no degrees of freedom left (James B. Grace et al. 2010). They represent a special class of model because they allow for everything to add up, meaning we can completely recover the observed matrix of covariances. Unsaturated models have testable implications, however. Under global estimation, our comparison for calculating the GOF is the saturated model. The reason is that saturated models permit every covariance to be explained, so the model fit function value \(F_{ML}\) fit function goes to zero. Thus, also \(\chi^2 = 0\). In comparison, an unsaturated model will have a positive \(F_{ML}\) (James B. Grace 2021).
Adding each a path from disk to richness
(rich) and evenness (even) turns our simple
model into a saturated one.
satur <-
"mass.above ~ nadd + rich + even + disk
rich ~ nadd + disk
even ~ nadd + disk
rich ~~ even"
fit.satur <- sem(satur, data = seabloom, estimator = "MLM")
summary(fit.satur, rsq = TRUE)## lavaan 0.6.17 ended normally after 32 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 12
##
## Number of observations 288
##
## Model Test User Model:
## Standard Scaled
## Test Statistic 0.000 0.000
## Degrees of freedom 0 0
##
## Parameter Estimates:
##
## Standard errors Robust.sem
## Information Expected
## Information saturated (h1) model Structured
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## mass.above ~
## nadd 0.096 0.007 14.018 0.000
## rich 0.106 0.019 5.639 0.000
## even -1.378 0.671 -2.054 0.040
## disk 0.446 0.073 6.079 0.000
## rich ~
## nadd -0.190 0.014 -13.573 0.000
## disk 0.412 0.269 1.534 0.125
## even ~
## nadd 0.006 0.001 8.503 0.000
## disk -0.003 0.009 -0.291 0.771
##
## Covariances:
## Estimate Std.Err z-value P(>|z|)
## .rich ~~
## .even -0.105 0.011 -9.614 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .mass.above 0.389 0.049 7.933 0.000
## .rich 5.201 0.514 10.122 0.000
## .even 0.006 0.001 10.693 0.000
##
## R-Square:
## Estimate
## mass.above 0.574
## rich 0.357
## even 0.297Exercise: What modification indices do you expect from a saturated model?
# modindices(fit.satur, minimum.value = 0)Model pruning
Now, let’s delete all statistically non-significant paths from the saturated model to obtain the most parsimonous model. There are two non-significant paths, i.e. between disturbance and richness and disturbance and evenness.
prune <-
"mass.above ~ nadd + rich + even + disk
rich ~ nadd
even ~ nadd
rich ~~ even"
fit.prune <- sem(prune, data = seabloom, estimator = "MLM")
summary(fit.prune, rsq = TRUE, fit.measures = TRUE)## lavaan 0.6.17 ended normally after 29 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 10
##
## Number of observations 288
##
## Model Test User Model:
## Standard Scaled
## Test Statistic 2.878 2.892
## Degrees of freedom 2 2
## P-value (Chi-square) 0.237 0.235
## Scaling correction factor 0.995
## Satorra-Bentler correction
##
## Model Test Baseline Model:
##
## Test statistic 592.561 765.397
## Degrees of freedom 9 9
## P-value 0.000 0.000
## Scaling correction factor 0.774
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.998 0.999
## Tucker-Lewis Index (TLI) 0.993 0.995
##
## Robust Comparative Fit Index (CFI) 0.998
## Robust Tucker-Lewis Index (TLI) 0.993
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -539.101 -539.101
## Loglikelihood unrestricted model (H1) -537.662 -537.662
##
## Akaike (AIC) 1098.202 1098.202
## Bayesian (BIC) 1134.831 1134.831
## Sample-size adjusted Bayesian (SABIC) 1103.120 1103.120
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.039 0.039
## 90 Percent confidence interval - lower 0.000 0.000
## 90 Percent confidence interval - upper 0.130 0.130
## P-value H_0: RMSEA <= 0.050 0.457 0.454
## P-value H_0: RMSEA >= 0.080 0.302 0.305
##
## Robust RMSEA 0.039
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.130
## P-value H_0: Robust RMSEA <= 0.050 0.456
## P-value H_0: Robust RMSEA >= 0.080 0.302
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.021 0.021
##
## Parameter Estimates:
##
## Standard errors Robust.sem
## Information Expected
## Information saturated (h1) model Structured
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## mass.above ~
## nadd 0.096 0.007 13.956 0.000
## rich 0.106 0.019 5.655 0.000
## even -1.378 0.665 -2.072 0.038
## disk 0.446 0.073 6.066 0.000
## rich ~
## nadd -0.190 0.014 -13.573 0.000
## even ~
## nadd 0.006 0.001 8.503 0.000
##
## Covariances:
## Estimate Std.Err z-value P(>|z|)
## .rich ~~
## .even -0.105 0.011 -9.600 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .mass.above 0.389 0.049 7.933 0.000
## .rich 5.243 0.536 9.786 0.000
## .even 0.006 0.001 10.754 0.000
##
## R-Square:
## Estimate
## mass.above 0.569
## rich 0.352
## even 0.297aictab(list(fit.satur, fit.prune),
c("saturated", "pruned"))##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## pruned 10 1099.00 0.00 0.68 0.68 -539.10
## saturated 12 1100.46 1.46 0.32 1.00 -537.66When comparing the CFI and AIC, we see that the pruned model is slightly superior to the saturated model. The AIC is the lowest (but with dAIC < 2), and the number of model parameters is lower and hence this model is the most parsimonous.
Summary output
Now that we have our final model, let’s have a closer look on the
model via summary requesting also the standardized
estimates:
summary(fit.prune, fit.measures = TRUE, rsq = TRUE, standardized = TRUE)## lavaan 0.6.17 ended normally after 29 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 10
##
## Number of observations 288
##
## Model Test User Model:
## Standard Scaled
## Test Statistic 2.878 2.892
## Degrees of freedom 2 2
## P-value (Chi-square) 0.237 0.235
## Scaling correction factor 0.995
## Satorra-Bentler correction
##
## Model Test Baseline Model:
##
## Test statistic 592.561 765.397
## Degrees of freedom 9 9
## P-value 0.000 0.000
## Scaling correction factor 0.774
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.998 0.999
## Tucker-Lewis Index (TLI) 0.993 0.995
##
## Robust Comparative Fit Index (CFI) 0.998
## Robust Tucker-Lewis Index (TLI) 0.993
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -539.101 -539.101
## Loglikelihood unrestricted model (H1) -537.662 -537.662
##
## Akaike (AIC) 1098.202 1098.202
## Bayesian (BIC) 1134.831 1134.831
## Sample-size adjusted Bayesian (SABIC) 1103.120 1103.120
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.039 0.039
## 90 Percent confidence interval - lower 0.000 0.000
## 90 Percent confidence interval - upper 0.130 0.130
## P-value H_0: RMSEA <= 0.050 0.457 0.454
## P-value H_0: RMSEA >= 0.080 0.302 0.305
##
## Robust RMSEA 0.039
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.130
## P-value H_0: Robust RMSEA <= 0.050 0.456
## P-value H_0: Robust RMSEA >= 0.080 0.302
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.021 0.021
##
## Parameter Estimates:
##
## Standard errors Robust.sem
## Information Expected
## Information saturated (h1) model Structured
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## mass.above ~
## nadd 0.096 0.007 13.956 0.000 0.096 0.897
## rich 0.106 0.019 5.655 0.000 0.106 0.316
## even -1.378 0.665 -2.072 0.038 -1.378 -0.137
## disk 0.446 0.073 6.066 0.000 0.446 0.235
## rich ~
## nadd -0.190 0.014 -13.573 0.000 -0.190 -0.593
## even ~
## nadd 0.006 0.001 8.503 0.000 0.006 0.545
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .rich ~~
## .even -0.105 0.011 -9.600 0.000 -0.105 -0.579
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .mass.above 0.389 0.049 7.933 0.000 0.389 0.431
## .rich 5.243 0.536 9.786 0.000 5.243 0.648
## .even 0.006 0.001 10.754 0.000 0.006 0.703
##
## R-Square:
## Estimate
## mass.above 0.569
## rich 0.352
## even 0.297Exercise: how would you interpret the meaning of the \(p\)-values of the single paths?
The meaning of the summary output is (James B. Grace 2021):
Degrees of freedomrepresents the number of paths omitted from the model, which provide a capacity to test the architecture of the modelp-value: probability of the data given our modelStd.err, the standard errorZ-value, the analogue to \(t\)-values derived from maximum likelihood estimation (ML)P(>|z|), the \(p\)-values, the probability of obtaining a \(z\) of the given value by chanceRegressions= the path coefficientsEstimates, the raw unstandardized coefficients
Variances= explained variance for endogenous variables = estimates of the error variancesEstimates, estimates of the error variances
R-squareis the variance explained for endogenous variables, the \(1 - error\) variance in standardized terms. Paths from error variables represent influences from un-modeled factors.
The inspect function retrieves information from a
lavaan object (for more options see the help):
inspect(fit.prune, what = "r2")## mass.above rich even
## 0.569 0.352 0.297Visualize the results
The package lavaanPlot
allows to simply and straight-forwardly visualize diagrams from
lavaan objects.
Another library that handles lavaan objects would be semPlot.
Then, we can plot the results with significance levels displayed as asterisks.
library("lavaanPlot")
lavaanPlot(model = fit.prune,
node_options = list(shape = "box", color = "gray",
fontname = "Helvetica"),
edge_options = list(color = "black"),
coefs = TRUE, covs = FALSE, stars = "regress")For more customization, check out the lavaanPlot
documentation.
# Add p-values and R^2: https://stackoverflow.com/questions/60706206/how-do-i-include-p-value-and-r-square-for-the-estimates-in-sempaths
# thePlot <- semPlotModel(fit.simple.up)
semPaths(fit.prune, what = "est", whatLabels = "est",
residuals = FALSE, intercepts = FALSE,
sizeMan = 10, sizeMan2 = 7, edge.label.cex = 1,
fade = FALSE, layout = "tree", style = "mx", nCharNodes = 0,
posCol = "#009e73ff", negCol = "#d55e00ff", edge.label.color = "black",
layoutSplit = TRUE, curve = 1, curvature = 1, #fixedStyle = 1,
exoCov = FALSE, rotation = 1)
What is missing in this plot are the \(R^2\) values for the endogenous variables. They can either be represented as the \(R^2\), the variance explained for endogenous variables or as the quantity of error variances (\(\zeta\)) either in raw or standardized units. Another option, if we wish to treat error variables like true causal influences, then we might use path coefficients for their effects. These are the square roots of the error variances (e.g., \(\sqrt{0.84} = 0.92\)) or alternatively, \(\sqrt{1 - R^2}\) (James B. Grace 2021).
Results to report
To present the results of a SEM, following should be stated in the report (James B. Grace 2021):
- Absolute GOF statistics for final model: \(\chi²\) with the affiliated \(p\)-value, CFI and degrees of freedom (df)
- Table of raw coefficients and statistics
- Table of total and indirect effects of interest
- Computed queries as table or graph